#include "CubicSpline.h"
#include <stdexcept>

// 构造函数
CubicSpline::CubicSpline(const QVector<double>& x, const QVector<double>& y) {
    if (x.size()  != y.size()  || x.size()  < 2) {
        throw std::invalid_argument("Invalid input data");
    }

    n = x.size()  - 1;
    this->x = x;
    this->y = y;
    computeSecondDerivatives();
}

// 插值函数
double CubicSpline::interpolate(double xi) {
    int i = findInterval(xi);
    double h = x[i + 1] - x[i];
    double a = (x[i + 1] - xi) / h;
    double b = (xi - x[i]) / h;
    return a * y[i] + b * y[i + 1] +
           ((a * a * a - a) * d2y[i] + (b * b * b - b) * d2y[i + 1]) * (h * h) / 6.0;
}

// 计算二阶导数
void CubicSpline::computeSecondDerivatives() {
    int size = x.size();
    QVector<double> h(size - 1), alpha(size - 1), l(size), mu(size), z(size);

    for (int i = 0; i < size - 1; ++i) {
        h[i] = x[i + 1] - x[i];
        if (i > 0) {
            alpha[i] = (3.0 / h[i]) * (y[i + 1] - y[i]) - (3.0 / h[i - 1]) * (y[i] - y[i - 1]);
        }
    }

    l[0] = 1.0;
    for (int i = 1; i < size - 1; ++i) {
        l[i] = 2.0 * (x[i + 1] - x[i - 1]) - h[i - 1] * mu[i - 1];
        mu[i] = h[i] / l[i];
        z[i] = (alpha[i] - h[i - 1] * z[i - 1]) / l[i];
    }

    l[size - 1] = 1.0;
    d2y.resize(size);
    d2y[size - 1] = 0.0;

    for (int i = size - 2; i >= 0; --i) {
        d2y[i] = z[i] - mu[i] * d2y[i + 1];
    }
}

// 找到 xi 所在的区间
int CubicSpline::findInterval(double xi) {
    if (xi <= x[0]) return 0;
    if (xi >= x[n]) return n - 1;
    for (int i = 0; i < n; ++i) {
        if (xi >= x[i] && xi <= x[i + 1]) {
            return i;
        }
    }
    return -1;
}
